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1  Motivation 

Initial  CFD  simulations  of  quasi-steady  hovering  flight  of  the  Drosophila 
sp.  have  studied  the  effects  of  delayed  rotation  of  the  wings  at  the  end 
of  each  stroke  and  it’s  effect  on  the  resulting  unsteady  aerodynamic  forces 
found  throughout  the  wing  [8].  Further  studies  were  carried  out  using  the 
same  CFD  model  as  before  to  determine  the  lift  and  power  requirements  for 
Drosophila  Virilis  [7].  To  verify  the  validity  of  the  simulation’s  results,  the 
specific  power  that  was  calculated  (the  sum  of  aerodynamic  and  intertial 
power  requirements  normalized  to  the  fruitfly  muscle  mass)  was  compared 
to  the  values  retrieved  by  Dickenson  on  tethered  fruitflies  [1].  Although  the 
periodic  behavior  of  the  ’translational’  and  rotational  motions  of  the  wing 
were  heavily  approximated,  the  comparison  showed  strong  agreement.  This 
lends  support  to  the  use  of  CFD  models  in  the  parametric  design  process 
used  in  the  pioneering  of  microscale  flapping  wings  [2].  It  is  of  great  im¬ 
portance  to  MEMs  designers  working  on  microscale  flapping  wing  to  study 
the  detachment  of  Leading  Edge  Vortices  (LEV).  Solving  the  Navier-Stokes 
equations  in  the  time- varying  coordinates  of  the  flapping  wing  proves  a  chal¬ 
lenge  computationally  but  offers  the  possibility  of  predicting  the  behavior 
of  the  LEV-shedding  phenomenon  to  which  flapping  wings  owe  a  majority 
of  their  lift.  In  addition  to  studying  the  shedding  of  LEVs,  CFD  simulations 
enable  one  to  quickly  assess  the  effect  of  wing  shaping  on  vertical  lift.  Pre- 
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liminary  code  was  developed  in  MATLAB  with  these  two  potential  studies 
in  mind. 


2  Governing  Equations 


The  code  that  was  developed  to  tackle  the  flapping  wing  problem  solves  the 
Navier-Stokes  equations  for  a  compressible  fluid  (written  here  in  cartesian 
coordinate  system): 
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where  u,  f,  g,  h,  fv,  gv,  hv 
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A  constant  density/ incompressible  fluid  assumption  can  be  applied  to 
equation  2  -  the  PDE  for  the  conservation  of  mass  -  which  sends  the 
term  to  zero.  The  resulting  PDE  can  also  be  written  as: 


V««  =  0  (10) 

This  statement  is  now  equivalent  to  saying  that  the  mass  flux  entering 
a  finite  volume  must  equal  the  mass  flux  leaving  that  volume.  The  reynolds 
number  for  the  flow  about  the  flapping  wings  of  Drosophila  Melanogaster 
was  determined  to  be  115  [3] .  The  low-reynolds  number  for  the  flow  allows  an 
additional  assumption  to  be  made:  the  flow  is  in  the  laminar  regieme.  This 
saves  us  from  the  added  complexity  of  implementing  a  model  for  turbulent 
fluid  flow.  The  following  form  of  the  Navier-Stokes  equations  is  developed 
by  Roger  and  Kwak  [4] 


t)  T 

[/tr+(^)n+1’m](£>n+1’m+1-£)n+1,m)  =  -Rn+1,m-^(1.5JDn+1’m-2Dn+0.5L»n-1) 

,!L>  '  (11) 

The  derivative  of  the  state  vector  with  respect  to  time  is  discretized  using 
a  second-order  three-point  backward  difference  formula.  The  distinguishing 
feature  of  this  equation  is  its  implementation  of  artificial  compressibilty. 

^  =  — /3V  •  un+1’m+1  (12) 

OT 

This  equation,  although  physically  meaningless  with  the  inclusion  of  the 
derivative  of  pressure  on  the  left  hand  side,  enforces  the  desired  incompress¬ 
ibility  assumption  which  lead  us  to  equation  10.  This  is  achieved  through 
carrying  the  form  of  the  Navier-Stokes  equations  forward  in  artificial  time. 

This  requires  adding  a  pseudotime  derivative  (a  derivative  with  respect  to 
artificial  time)  to  all  four  of  the  Navier-Stokes  equations.  The  three  other 
equations  that  complement  equation  12  can  be  written  in  a  similar  fashion 
to  the  one  reproduced  below  (in  its  discretized  form). 


1.5 un+hm+l  _  L5fin+l,m  _  _p+l,m+l  _  1.5 +  ^ 

At  =  At  ^13) 

Together,  equation  12  and  the  three  equations  that  take  on  the  form 
of  13  can  be  written  compactly  in  the  delta  form  (equation  1)  mentioned 
earlier: 

where  D,  R,  E,  F,  G,  Ev ,  Fv,  Gv  are 
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and  where  dR/dD  consists  of  the  jacobians  of  the  convective  and  viscous 
flux  terms.  The  jacobian  arises  from  a  common  technique  used  in  solving  the 
Navier-Stokes  equations  numerically.  This  technique  linearizes  the  residuals 
of  each  equation  so  that  one  can  write  the  equations  in  implicit  form  and 
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from  this  the  systems  of  linear  equations  can  be  solved  numerically  using 
conventional  methods. 

3  Numerical  Method 

Equation  11  is  written  in  a  form  that  can  be  solved  using  an  procedure 
called  the  Gauss-Seidel  line-iteration.  The  viscous  fluxes  are  discretized 
with  a  second  order  central  difference  scheme  while  the  convective  fluxes 
are  discretized  using  a  cell-centered  scheme  that  treats  the  fluxes  across 
each  node  (each  finite  volume)  as  conservative  variables.  This  discretization 
proves  much  more  stable  as  it  provides  numerical  dissipation  automatically. 
This  discretization  is  formed  as  follows  for  the  flux  passing  through  the  faces 
corresponding  to  x-ward  normal  vectors: 

[/t +1/2  ~  ft- Vg]  _  q 

Ax 

fi+ 1/2  =  \[f(Qi+ 1)  +  /(?<)]  -  \[^f?+ 1/2  -  A/"+l/2]  (23) 

A/?+i/2  =  A±(q)Aqi+ 1/2  (24) 

the  last  equation  24,  is  the  flux  difference  term  found  in  equation  23 
and  it  can  be  calculated  provided  that  one  can  find  a  matrix  A  that  is 
the  jacobian  of  the  flux  vector  with  respect  to  the  state  vector  D.  This 
jacobian  must  also  carry  with  it  other  properties  all  captured  by  the  so 
called  ’property  U’.  Creating  a  jacobian  matrix  with  ’property  U’  allows  the 
jacobian  to  be  representative  of  the  change  in  flux  with  respect  to  state 
vectors  at  the  interface  between  the  two  nodes  it  is  responsible  for.  P.  L. 
Roe  developed  a  method  to  produce  such  jacobians  that  is  implemented  on 
the  Euler  equations  [3]  but  can  be  partly  extended  toward  the  navier-stokes 
equations  and  its  convective  fluxes. 

4  Results  and  Further  Work 

The  main  problem  facing  the  code  now  in  its  development  is  that  upon 
attempting  to  solve  the  linear  system  of  equations  according  to  the  Gauss- 
Seidel  line-iteration,  the  matrix  cannot  be  numerically  inverted  because  it  is 
found  by  linsolveQ  in  MATLAB  to  be  singular.  This  prevents  results  from 
being  produced  for  the  simple  3D  channel  problem. 


(22) 
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